#######################################################################################
# Summary of Simulation Results: Monte Carlo Study 3 
#################################################################################

# Function to Average Across S=1000 Simulations

MeanOutFun <- function(IN_FILE, OUT_FILE){

result.array <- readRDS(paste(IN_FILE, ".rds", sep=""))
mean.array <- matrix(NA, 6, 8)
rownames(mean.array) <- rownames(result.array)
colnames(mean.array) <- c("NumRob", "TP", "FP", "TN", "FN", "CSign", "RMSE", "Corrs")
for(i in 1:6){
  for(j in 1:8){
    
    mean.array[i, j] <- mean(result.array[i, j, ], na.rm=T) 
    
  }
}

#write.table(round(mean.array, 2), file=paste("/Users/richardtraunmuller/Dropbox/PSRM_Replication/", OUT_FILE,  ".csv", sep=""), sep=",")
print(mean.array)
}

##################################################################################################################################
# Table 3: Probability of False Positive and False Negatives for Varying Levels of Correlation 

p <- 5

n1 <- round(1-MeanOutFun("five_norm_1_result_mat_1", "five_norm_1_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("five_norm_1_result_mat_2", "five_norm_1_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("five_norm_1_result_mat_3", "five_norm_1_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("five_norm_1_result_mat_4", "five_norm_1_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("five_norm_1_result_mat_5", "five_norm_1_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("five_norm_1_result_mat_1", "five_norm_1_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("five_norm_1_result_mat_2", "five_norm_1_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("five_norm_1_result_mat_3", "five_norm_1_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("five_norm_1_result_mat_4", "five_norm_1_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("five_norm_1_result_mat_5", "five_norm_1_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.3.a <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

#
p <- 5

n1 <- round(1-MeanOutFun("five_norm_3_result_mat_1", "five_norm_3_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("five_norm_3_result_mat_2", "five_norm_3_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("five_norm_3_result_mat_3", "five_norm_3_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("five_norm_3_result_mat_4", "five_norm_3_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("five_norm_3_result_mat_5", "five_norm_3_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("five_norm_3_result_mat_1", "five_norm_3_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("five_norm_3_result_mat_2", "five_norm_3_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("five_norm_3_result_mat_3", "five_norm_3_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("five_norm_3_result_mat_4", "five_norm_3_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("five_norm_3_result_mat_5", "five_norm_3_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.3.b <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

#
p <- 5

n1 <- round(1-MeanOutFun("five_norm_5_result_mat_1", "five_norm_5_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("five_norm_5_result_mat_2", "five_norm_5_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("five_norm_5_result_mat_3", "five_norm_5_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("five_norm_5_result_mat_4", "five_norm_5_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("five_norm_5_result_mat_5", "five_norm_5_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("five_norm_5_result_mat_1", "five_norm_5_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("five_norm_5_result_mat_2", "five_norm_5_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("five_norm_5_result_mat_3", "five_norm_5_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("five_norm_5_result_mat_4", "five_norm_5_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("five_norm_5_result_mat_5", "five_norm_5_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.3.c <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

table.3 <- rbind(tab.3.a, tab.3.b, tab.3.c)
rownames(table.3) <- rep(rep(c("EBA_Leamer", "EBA_sign", "EBA_90", "MA_unig", "MA_Sala", "MA_BIC"), each=2), 3)
colnames(table.3) <- c(".25", ".50", ".75", "1.00", "1.25")

write.table(table.3, file="pluemper_traunmueller_table3.csv", sep=",")

print(table.3)
############################################################################################################################


